Time-Dependent Gutzwiller Theory for Multiband Hubbard Models 
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Based on the variational Gutzwiller theory, we present a method for the computation of response 
functions for multiband Hubbard models with general local Coulomb interactions. The improvement 
over the conventional random-phase approximation is exemplified for an infinite-dimensional two- 
band Hubbard model where the incorporation of the local multiplet-structure leads to a much 
larger sensitivity of ferromagnetism on the Hund coupling. Our method can be implemented into 
LDA+ Gutzwiller schemes and will therefore be an important tool for the computation of response 
functions for strongly correlated materials. 
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During the last two decades a number of schemes has 
been developed to improve density functional related the- 
ories (DFT) [l|, [2| for materials with medium to strong 
Coulomb-interaction effects. For such systems, the stan- 
dard DFT methods are usually insufficient since they are 
based on effective single-particle approaches. In this con- 
text, a common idea is to separate the electrons in local- 
ized orbitals (typically those in d- or f-bands) from the 
delocalized ones and to assume for the former a Hubbard- 
type interaction which is then treated in a given approx- 
imation. In the LDA+U method [3| the local-density 
approximation (LDA) is combined with a Hartree-Fock 
(HF) like treatment of the local interaction, whereas in 
the LDA+DMFT approach [4| dynamical mean-field the- 
ory (DMFT) is applied to deal with the Hubbard term. 
The LDA+DMFT method is quite successful in calculat- 
ing single-particle excitations of correlated materials [3[ ; 
however, it is rather demanding from a numerical point 
of view since it requires the self-consistent solution of 
complex single-impurity models. In this regard, the 
Gutzwiller variational theory is a promising alternative 
to be combined with LDA computations since it joins the 
numerical simplicity of an LDA+U calculation with the 
accuracy of DMFT methods. The Gutzwiller theory is 
based on an exact evaluation of Gutzwiller wave func- 
tions in the limit of infinite spatial dimensions D. This 
has been accomplished for general multiband Hubbard 
models in Rcf. [6(. Applying the infinite D energy func- 
tional to finite dimensional systems is usually denoted 
as the 'Gutzwiller approximation' (GA). The multi-band 
GA has been successfully combined with bandstructure 
calculations and, among others, applied to ferromag- 
netic Nickel @, M, and LaOFeAsNa Q. Self-consistent 
Gutzwiller-DFT schemes have been developed and ap- 
plied in Refs. [lMi- 

A comparison of experimental data with results from 
DFT often requires the computation of response func- 
tions where one is interested in the frequency and wave- 
vector dependence of two-particle correlation functions. 
Within DFT, such excitations can be computed via 
the time-dependent DFT (TDDFT), which is based on 



a proper generalization of the Hohenberg-Kohn theo- 
rem [15j. In general, when the external perturbation 
is small the susceptibility XgO^) °f the interacting elec- 
trons can be expressed via the fictitious response function 
of the Kohn-Sham reference system. The structure of 
the resulting self-consistency equations for x g (w) bears 
strong resemblance with the random-phase approxima- 
tion (RPA) where the interaction kernel is now composed 
of the Hartree- and exchange-correlation potential. 

The difficulties of the DFT to describe materials 
with strong Coulomb interactions equally concern the 
TDDFT. Therefore, the question arises whether the con- 
cept of TDDFT can be applied to the aforementioned 
improvements of DFT for such systems. Concerning 
DMFT, there are obvious limitations due to the numer- 
ical complexity of this approach. On the other hand, 
it has been shown that a time-dependent GA (TDGA) 
can be constructed [16| for single-band Hubbard models. 
This approach has been successfully applied to a num- 
ber of systems and response functions [17H20J . In order 
to construct a time-dependent theory which can be com- 
bined with Gutzwiller-DFT computations, the concept of 
the TDGA has to be generalized to multiband Hubbard 
models with arbitrary local interactions. This is the main 
purpose of the present paper. 

We consider a general class of multi-band Hubbard 
models which is defined as 
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(1) 



Here, the first term describes the hopping of electrons 
between N spin-orbital states a, a' on L lattice sites i,j, 
respectively. The Hamiltonian H\ oc _i contains all local 
terms, i.e., two-particle Coulomb interactions and or- 
bital onsite-energies. We further introduce the eigen- 
states \Yi) of H\ OC i and the corresponding energies E^"?, 

-ffloc,i|rj) = E^jlTi) . 

Within the Gutzwiller theory [5|, |6|, the Hamilto- 
nian (JlJ is investigated by means of the variational wave 
function |* G ) = P G |*o) = Eh-Pil^o) , where |* > is 
a normalized single-particle product state and the local 



i.e 



Gutzwillcr correlator Pj = 2 r r » Ap r , |r),-,-(r'| depends 

on the matrix A of variational parameters Ar,r' • 

The expectation value of the Hamiltonian ([1]) in infinite 
dimensions (i.e., in GA) is read as 



£ GA (p,A) = Yl 4"^ V ^J Pi*>,io 

+ L J2 ^f° CA f,r A r,f' m f,f' ( 2 ) 
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where 



((|r)(r'|))$ can be expressed [6( via the 
local part of the uncorrelated density matrix p with the 



i T _ T , — 



elements 



Pjc 



1 \Ci,<J C j,a')y< 



The 'renormalization 



factors' q°, in ^ account for the correlation induced 
band narrowing and depend on A and p. 

The energy functional Eq. ([5]) has to be minimized 
with respect to the n var variational parameters Ar.r' 
and the elements of the density matrix p [2l|. Within 
the framework of the GA the variational parameters 
have to obey the completeness condition C\ (p,X) = 
(Pi Pj)* — 1 = plus n con further local constraints 

C^ a/ (p,X) = Pij^ - (P % A4<V)*o = which Sug- 
gests the definition of the Lagrange function 

£ GA (p,A) =£ GA (p,A)+]>>r^%A) 

+ E A ^^(P,A). (3) 

i : (7(7' 

The GA ground state is thus defined from the condition 
SC = which fixes all n var variational parameters Ar.r' 
plus the ri con + l Lagrange parameters A = {A^ , A^ j CT ,}. 
Having thus determined the GA ground state, one can de- 
fine an auxiliary 'Gutzwillcr single-particle Hamiltonian' 
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(4) 



the diagonalization of which yields also the unoccupied 
single-particle orbitals in addition to the occupied ones 
already implicit in the ground-state density matrix p°. 

A small time-dependent perturbation will induce small 
(time-dependent) deviations in both the density matrix 
(dp) and in the variational and Lagrange parameters (SX, 
S A) from their ground state values. Following the spirit of 
time-dependent HF theory (TDHF) [see e.g. [13, Q] the 
response of the system can be obtained from an expansion 
of the energy functional up to second order in the density 
fluctuations. The complication in the present case is due 
to the additional fluctuations in the variational parame- 
ters whose time-dependence is unknown. In case of the 
one-band model one can use the constraints to express 
all variational parameter fluctuations but one (which is 
usually chosen as the double occupancy fluctuation Sd) 
via the density fluctuations 5p (and Sd). In the multi- 
band case studied here, one usually has n var 3> n con 



and there is no general analytical way to express some 
"■con + 1 dependent variational parameters via the re- 



maining n, 



1 independent ones and the density 



matrix. This problem can be circumvented by the fact 
that the expansion of E G in the independent parameter 
fluctuations is equivalent to an expansion of the Lagrange 
function Eq. ([3]) in all the fluctuations Sp and SX. Ex- 
panding Eq. ([3]) up to second order yields 
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where h° is given by the GA ground-state result for (j4j). 
For clarity, we have casted the multiplicity of dependen- 
cies of the fluctuations in a single index. Like in the 
single-band case [la ] we can now relate the dynamics 
of the n var + n con + 1 parameter fluctuations SX and 
(5A to that of the density fluctuations Sp requiring that 



jr- = and "~. - = 0. The resulting linear system 

aoX ooA ^ 

of equations can be inverted yielding SX = 5X[Sp] and 
<5 A = SA[Sp] which upon inserting into Eq. (j3|) results 
in a Lagrange function which depends on density matrix 
fluctuations only SC GA (Sp) = Tr(h°Sp) + \K,^SpiSpj . A 
more general approach which explicitely incorporates the 
dynamics of the variational parameters may be derived 
along the lines presented in [24}. We expect, however, 
that this improvement will change the results mostly at 
high energies. 

The remaining steps are analogous to TDHF [22J, |23 | 
and TDGA for the one-band model [l6|. We construct 
the 'bare' susceptibilities X?.-(w) on the GA level from 
the eigenstates of the Hamiltonian h°. The correspond- 
ing correlation functions of the interacting system are 
then obtained from the RPA-like series Xij{ u ) = X ( !j( UJ ) + 
X° n (w)/C nm Xmj(w). Note that , in contrast to conven- 
tional TDHF theory the TDGA induces also couplings 
between transitive density fluctuations (i.e., ~ 3( c i a c j<?')) 
which enlarges the size of the RPA-like matrix equation. 
Nevertheless, the numerical efforts of our method ex- 
ceed those of an ordinary HF+RPA only marginally. It 
is also straightforward to implement our approach into 
LDA+Gutzwiller schemes outlined in Refs. 0, [tHH. 
The number of parameters Xi needed in the expansion 
((5|) depends on the symmetries of the model and the re- 
sponse function one aims to calculate. In our calculations 
of magnetic susceptibilities for a two-band model (see be- 
low) 16 such parameters have been used. 




FIG. 1: (color online) Frequency and momentum dependence 
of the magnetic susceptibility for finite ferromagnetic ground- 
state magnetization n-|-— nj. = 0.677 per orbital, various values 
of rj, doping 8 = 0.6, and U/t = 10.5, J/U = 0.2. The main 
panel covers the frequency range up to the Stoner excitations 
whereas the inset shows the low energy magnon part. 



We exemplify our method for a Hamiltonian with two 
degenerate bands on a hypercubic (he) lattice (intra- 
orbital hopping tij = t/y/2T>) where the GA is ex- 
act for the Gutzwiller variational wave-function. Note 
that for the he lattice the momentum dependence of 
the susceptibilities is contained in a single par ameter 
V = ?7q = ^[cos(gi) + cos(q 2 ) + .... + cos(<7d)] [26| which 
takes values — 1 < r\ < +1. We work with the standard 
atomic Hamiltonian for two degenerate e g orbitals [6j 
with intra-orbital Coulomb interaction U and exchange 
interaction J. Note that in contrast to most previous 
DMFT investigations of similar models (cf. Ref. [25| and 
references therein) we explicitly include the pair hopping 
term to guarantee rotational invariance in orbital-spin 
space. Fig. Q] displays the spin susceptibility x(q, w ) = 
f dtexp(icjt)(TSq(t)S~(0)) for a ferromagnetic system 
(magnetization m = 0.677 per orbital) where S+ denotes 
the spin-flip operator in momentum space. The spectrum 
is composed of a low energy magnon (cf. inset to Fig. [1]) 
and the higher energy Stoner contribution. In the limit 
q — > (rj — > 1) and u) — > one obtains the Goldstonc 
mode as a delta-function (broadened by e = 10 _5 £) in 
Fig.[I} which comprises all the spectral weight and which 
provides an additional consistency check of our approach. 
For finite momenta the magnon excitations get rapidly 
damped (cf. inset to Fig. [T]) due to the large phase space 
of spin flip particle-hole decay processes in infinite di- 
mensions. Contrary to HF+RPA, where the Stoner exci- 
tations are expected to appear at u ~ (U + J)m (= 8.53i 
for the same magnetization) one observes that these are 
reduced to ~ 1.22i within the TDGA thus rapidly merg- 
ing with the magnon excitations for rj < 1. With regard 
to this correlation induced reduction of the Stoner pa- 
rameter, our findings are similar to local Ansatz inves- 
tigations |27j but they generalise them by providing us 
with the full spin-excitation spectrum. 

Figs. [5] and [3] report the magnetic phase diagrams for 
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FIG. 2: (color online) TDGA magnetic phase diagram for the 
two-band model on the hypercubic lattice. Shown are the 
phase boundaries between paramagnetic (PM), incommensu- 
rate (IC) and stable ferromagnetic (FM) phases. The full 
square at 5 = 1 indicates the BR transition towards a local- 
ized regime (thick solid line). The solid circle marks a Lifshitz 
point where PM, IC, and FM phases merge. AF order is sta- 
ble at 5 = (solid green line). In panels (a,c) also the orbital 
ordered (OO) phase is shown (dashed dotted). Doping 5 is 
defined as the concentration of holes per site with 5 = 0(1) 
being the state with 2(1) particles per site. 



both TDGA and HF+RPA and various values of J/U, 
respectively. We have determined the second order tran- 
sition between paramagnetic (PM) and incommensurate 
(IC) or PM and ferromagnetic (FM) phases from the con- 
dition DET|1 - X a in (u = 0)/C„ m | = 0. For the PM-FM 
transition we verified the second-order nature by calcu- 
lating the ground-state energy as a function of the mag- 
netic moment. Especially the evaluation of IC phase 
boundaries is straightforward in our approach whereas 
the local nature of DMFT has hampered the determina- 
tion of these phases in related investigations 28l43fj| . The 
solid lines in Figs. [2] and [3] are the enveloping curve for 
all the phase boundaries with different 7y q . At half filling 
(S = 0) an antiferromagnetic (AF) phase is stable over 
the whole range of U/t (rj = —1) due to perfect nest- 
ing. It evolves continuously into IC ( — 1 < ?/ < 1) and 
eventually FM (rj = 1) phases at finite doping. Stabil- 
ity of the latter has been determined from the magnetic 
excitation spectrum (cf. Fig. [1]). The dashed lines in 
Figs. [2] and [3] separate a stable FM state with positive 
'spin-wave stiffness' (i.e. a positive slope of the magnon 
dispersion in the ferromagnetic state) from the IC phases 
where a FM state would decay due to quantum fluctu- 
ations. Note that due to the stability of the AF phase 
at S = these phase boundaries diverge as a function of 
U/t upon approaching half-filling. 

The HF+RPA phase diagram shows the known limi- 
tations of a mean-field multiband approach: It depends 



only weakly on J/t and it contains a stable FM state in 
the whole doping regime at sufficiently large U/t. In addi- 
tion, one always finds a a Quantum Lifshitz point (QLP) 
where IC, PM, and FM phases meet, similar to previous 
investi gati ons in the one-band Hubbard model on the he 
lattice [23| . In contrast, within the TDGA a stable FM 
phase requires a minimum value of J/U (cf. Fig. [2ji) 
which favors local triplet formation in agreement with 
DMFT J2J|. In the intermediate J/U regime (Figs. Eb,c) 
IC and FM phases are bounded by the Brinkman-Ricc 
(BR) localization transition which occurs at 5 = 1 and 
U c /t (the latter depending on J/U). Note that the same 
termination of the FM phase has been found in DMFT in- 
vestigations on a three-orbital model [30( • Mapping to a 
Kugcl-Khomskii Hamiltonian [3l[ suggests the formation 
of additional orbital order (00). This 00 is expected 
around doping 8 = 1 being FM in the spin- and AF in 
the orbital channel, i.e. (o f), (t °), (° t), (t °)- DMFT 
investigations of this issue within the two-orbital model 
remain ambiguous [28|, |29(- In Fig. |2K, c the 00 phase 
as obtained from minimizing the GA functional is shown 
which in the PM regime occurs as a first order transi- 
tion. Within the FM phase the transition towards 00 
is second order for J/U = 0.2. However, we find that it 
can also be of first order for larger values J/U = 0.3 (not 
shown). For all values of J/U the BR transition at S = 1 
is masked by 00. It thus appears as a generic finding that 
PM insulating phases at integer doping are prevented by 
another instability, either magnetic (AF at S = 0) or due 
to 00 (at 5 = 1), similar to the three-orbital model |30|. 

The multiband TDGA can be easily generalized to the 
pair and charge sector where in the latter it allows for the 
computation of response functions like optical conductiv- 
ity and polarizability for multiorbital Hubbard models. 
It is obvious that the formalism offers a valuable tool with 
manageable numerical effort for the investigation of the 
dynamics of multiband Hubbard models in the context 
of correlated materials. 
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FIG. 3: Same as Fig. |5] but obtained within HF+RPA. 



